function y=bayesmixture2(Data,Demo,mixedlogitbeta,K)

global H

S=10000;                       %number of simulation

N=size(Data,1)/40;                              %number of agent

logitbeta=mixedlogitbeta;
J=(size(logitbeta,1)+1)*size(Demo,2);           %number of parameters
D=size(Demo,2);                                 %number of demo var
option=max(Data(:,3));
L=size(Data,2)-4+option-1;

%Prior
rhoR=(2.93/(sqrt(L)))^2;
rho=(2.93/(sqrt(L*D)))^2;
a=.4*N/K;                           %prior on alphas in Dirichlet
b_bar=zeros(K,L);
v=L+3;
V=v*eye(L,L);                             %prior on sigma
A_beta=0.2*eye(D*L);
A_b=.1;
beta_bar=zeros(D*L,1);


%Set intial values for the simulation
%pvec=initialpvec;
pvec=ones(1,K)./K;
rbeta=mvnrnd(zeros(N,L),2*eye(L));
beta=zeros(D,L);
Beta=Demo*beta+rbeta;
ind=randsample(K,N,true,ones(1,K)/K);

for k=1:K
    W(:,:,k)=3*eye(L);
    b(k,:)=zeros(1,L);
end


for s=1:S
    %Start Gibbs Sampling
    %Step 1: Draw ind | pvec, mu_k, W_k, Beta
    Beta_star=Beta-Demo*beta;
    %Set up the prob in the multinomial for each n
    for k=1:K
        p(:,k)=pvec(k)*mvnpdf(Beta_star,b(k,:),W(:,:,k));
        p(p(:,k)==0,k)=1e-299;
    end
    
    denom=sum(p,2);
    pi=p./(denom*ones(1,K));
    for n=1:N
       ind(n)=randsample(K,1,true,pi(n,:));
    end
    clear denom pi p

    %Step 2: Draw pvec | ind
    for k=1:K
        Kvec(1,k)=k;
    end

    ncomp=histc(ind,Kvec)';
    alpha=ncomp+a;

    %draw pvec from Dirichlet
    for k=1:K
        p(k,:)=gamrnd(alpha(k),1);
    end
    for k=1:K
        pvec(k)=p(k)/sum(p);
    end
    clear p
    
    %Step 3: Draw W and b using multireg

    for k=1:K
       temp=multireg(Beta_star(ind==k,:),ones(ncomp(k),1),v,V,b_bar(k,:),A_b);
       b(k,:)=temp(1,:);
       W(:,:,k)=temp(2:size(temp,1),:);
       clear temp
    end
    
    %Step 4: Draw beta
    
    %Stack and compute some moments needed for simulation
    XX=kron(zeros(D,D),V);
    Xy=inv(V)*zeros(N,L)'*zeros(N,D);
    for k=1:K
       XX=XX+kron(Demo(ind==k,:)'*Demo(ind==k,:),inv(W(:,:,k)));
       Xy=Xy+inv(W(:,:,k))*(Beta(ind==k,:)-ones(ncomp(k),1)*b(k,:))'*Demo(ind==k,:);
    end
    Xy=reshape(Xy,[],1);
    beta=inv(XX+A_beta)*(Xy+A_beta*beta_bar);
    beta=(reshape(beta,L,D))';
    
    %Step 5: Use RW MH algorithm to draw Beta
    count=0;
    for n=1:N
        Wn=rho*inv(H(:,:,n)+inv(W(:,:,ind(n))));
        
        Datan=Data((n-1)*40+1:n*40,:) ;
        
        %Draw Beta1 from multivariate normal with mean Beta(n,:) and
        %variance Wn
        Beta1=Beta(n,:)+(chol(Wn)'*randn(L,1))';
        
        %ranklikelihood([rbeta(n,:) beta1], n, Data, idcharac)
        posterior=(likelihood(Beta(n,:),Datan)*mvnpdf(Beta(n,:)-Demo(n,:)*beta,b(ind(n),:),W(:,:,ind(n))));
        if posterior==0
            Beta(n,:)=logitbeta(1:L);
            disp('reset initial values')
        elseif isnan(posterior)==1
            Beta(n,:)=logitbeta(1:L);
            disp('reset initial values')
        else
            mu=rand;
            F=likelihood(Beta1,Datan)*mvnpdf(Beta1-Demo(n,:)*beta,b(ind(n),:),W(:,:,ind(n)))/posterior;
            if mu<=F
                Beta(n,:)=Beta1;
                count=count+1;
            end
        end
        clear Datan
    end

    %Display s, ncomp, mean(Beta) while running the program
    s
    count
    ncomp
    beta
    mean(Beta)
    %save the estimates of Betas, W(:,:,k), b(k,:), pvec 
    y(s,:)=[ncomp ind' pvec reshape(Beta',1,[]) reshape(beta,1,[])...
         reshape(b',1,[]) reshape(W,1,[])];

end
% 
% 
